Superconductivity in graphene stacks: from the bilayer to graphite 
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We study the superconducting phase transition, both in a graphene bilayer and in graphite. For 
that purpose we derive the mean-field effective potential for a stack of graphene layers presenting 
hopping between adjacent sheets. For describing superconductivity, we assume there is an on-site 
attractive interaction between electrons and determine the superconducting critical temperature as 
a function of the chemical potential. This dis plays a dome-shaped curve, in agreement with previ- 
ous results for two-dimensional Dirac fermion a ^1 We show that the hopping between adjacent 
layers increases the critical temperature for small values of the chemical potential. Finally, we con- 
sider a minimal model for graphit^^ and show that the transition temperature is higher than that 
for the graphene bilayer for small values of chemical potential. This might explain why intrinsic 
superconductivity is observed in graphite. 



PACS numbers: 74.25.Dw,74.70.Wz 

I. INTRODUCTION 

Graphene is a one-atom-thick layer of graphitd^. The 
carbon atoms in each layer are arranged in a honey- 
comb lattice and the tight-binding energy presents a band 
structure such that the valence and conduction bands 
touch precisely in the vertices of two inequivalent Dirac 
cones in the Brillouin zone. The electronic excitations 
appearing in the conduction band have the dispersion 
relation of a relativistic massless particle and their prop- 
erties, accordingly, will be determined by the Dirac equa- 
tion. Graphene is believed to be the parent compound of 
most of the carbon-based systems and their electric, mag- 
netic and elastic properties all originate from the prop- 
erties of graphene. 

Interestingly, several carbon-based compounds present 
superconductivity. For instance, the graphite interca- 
lated compounds (GlCj^ which consists of graphene 
sheets alternated by alkali layers, mainly acting as charge 
reservoirs, becomes superconducting with the transition 
temperature ranging from below IK for KCg to 11.5 K 
for CaC 6^"^; some fuUerides present critical tempera- 
tures as high as 33 K as applied pressure or the chemical 
composition increases the lattice parametei^^ and there 
are reports of room temperature local superconductiv- 
ity within isolated "grains" in highly oriented pyrolitic 
graphite (HOPGj^ and also with critical temperature 
Tc ~ 25 K in thin sampled. Moreover, a fully saturated 
hydrocarbon derived from a single graphene sheet, called 
graphane, is predicted to be a high-temperature electron- 
phonon superconductor exhibiting a critical temperature 
of above 90 K^. 

Despite the fact that theoretical conjectures have been 
proposed as possible candidates to produce superconduc- 
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intrinsic superconductivity has never been ob- 
served in graphene, but it could be only induced by prox- 
imity effects, where a superconducting current propa- 
gated through a superconductor-normal-superconductor 
(SNS) Josephson junction, with graphene as the N re- 
gion^®. 

Nevertheless, the stability of the superconducting 
phase has been investigated in graphen^i^''^ and the 
symmetry of the order parameter in the honeycomb lat- 
tice was identified; if there is an on-site net attractive 
interaction between electrons in the honeycomb lattice, 
the usual s-wave singlet pairing is favourecP^l. As nearest- 
neighbours attraction are taken into account, an exotic 
combination of s-wave and p-wave superconducting or- 
der parameters is possibl^^H jn the context of the t-J-U 
model, /-wave triplet-pairing and d -\- id singlet-pairing 
instabilities are found to emerge away from half-filling^J^. 

Previously, some of us have investigated the phase di- 
agram of a quasi-two-dimensional interacting Dirac elec- 
trons system forming Cooper pairs in the singlet state, 
which is a suitable model to describe a stack of uncoupled 
superconducting graphene sheets, and we have found a 
quantum critical point connecting the normal and super- 
conducting phases at a certain critical couplinj^l. If low 
magnetic fields are applied to the system, we have found 
a critical field as a function of the superconducting inter- 
actioiP^. 

In those previous investigations, the variation of the 
chemical potential was not taken into account; however, 
applying a bias voltage, the carrier density of graphene 
can be controlled by electric field effect. Therefore, in 
the present paper we investigate the effect of the chem- 
ical potential as a free parameter of our model and we 
also consider the effect of the out-of-plane hopping be- 
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tween adjacent graphene sheets. In order to describe 
graphite, we consider the minimal model with the elec- 
tron tunneling between the nearest sites in the plane and 
out of the plane. We have found that the superconduct- 
ing critical temperature is enhanced at small values of 
the chemical potential for graphite when compared to 
the values predicted by us for graphene bilayer, what 
might explain why intrinsic superconductivity has been 
observed in HOPG. 

The paper is organized as follows: In Sec.|ll]we present 
the model Hamiltonian for the graphene bilayer, the dis- 
persion relation is calculated and the effective potential 
(free energy) is derived. In Sec. II A the superconducting 
phase diagram at T = is obtained analyzing the min- 
ima conditions for the effective potential for several val- 
ues of the interaction and the hopping between adjacent 
graphene sheets. In Sec. II B we calculate the supercon- 
ducting critical temperature as a function of the chemical 
potential for several values of the hopping parameter be- 
tween layers. The results represent an upper bound for 
the Kosterlitz-Thouless transition. In Sec. lIIII our results 
for the superconducting phase diagram are extended for 
an infinite number of coupled graphene layers consider- 
ing the electron tunneling amplitudes between the near- 
est sites in the plane and out of the plane, 
devoted to the conclusion. 



II. GRAPHENE BILAYER 

Consider a stack of N graphene layers with a hopping 
term between adjacent planes, where the upper layer has 
its B sublattice on top of sublattice A of the underlying 
layer (Bernal stacking), as can be seen in Fig. [l] The 
Hamiltonian of each coupled layer is described by the 
following-^3 
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where the index I = 1, - ■ ■ ,N characterizes the different 
planes and ji is the chemical potential. The second line in 
the RHS of the above equation describes the hopping be- 
tween electrons of different sublattices within a graphene 
sheet, while the third line describes the hopping between 
layers. The hopping parameter is about t « 2.8 eV and 
t_L ~ t/10. The operators a\„i = Ylik^^^^^ a i ^'^'^ 
b\^i = X)k e-"^^'^' fejc ; create, respectively, an electron 
on site i with spin a on sublattice A and an electron on 
site i with spin a on sublattice B of plane I. In the hon- 
eycomb lattice we have Sfc = 1 + e^^ '^^ + e^^'^'^ , where 
ai — acx and 2a.2 — a (e^, — \/3ej,) , as shown in Fig.jl] 
The lattice parameter is a = 2.46 Afor graphene. 



^ R9 




Figure 1: Lattice structure of two adjacent graphene layers 
(aftePl). 



We add an on-site attractive interaction between the 
electrons within each graphene layer forming Cooper 
pairs in the s-wave state. The interaction term is given 
by 



sc,i 
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— k', — cr,/'^k' cr.j 
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k,CT,r-k -<7,i^-k' -(T,i&k',<T,ij , (2) 

with g > 0. The origin of the interaction is to be de- 
termined by some underlying microscopic theory, which 
is not considered here. However, the symmetry of the 
gap originated from this interaction is consistent with the 
isotropic s-wave symmetry gap observed in some GIC^^Zl, 
Introducing the following Nambu fermion field. 



where 



one can rewrite the combined Hamiltonian Ht^i + Hsc,i 
at the mean-field level. 
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where, by definition, the superconducting order parame- 
ter is 



^ k k 

and the 8x8 matrix A in Eq. ([5| is given by 
A 

with 



Ai Ai2 
A21 A2 



(6) 



(7) 



Ai ^A2 



I —[Jl —tsk A ' 

-tsl -fi A 

A* n tsl 

\ A* tsk II 



(8) 



3 



and 
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From i?MF in Eq. ([5|, follows the dispersion relation, 
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Figure 2: Band structure for (a) A = /i = 0, (b) A = and 
= 0.2t, (c) A = O.lt and ^ = (d) A = O.U and /i = 0.2f. 
Energy is given in units of t and = 0.2t. 

Let us neglect the hopping term between planes for a 
moment and consider only the normal state of the sys- 
tem at the Fermi level, which means — IS. = \x = 
In that case, A has eight eigenvalues, but only two are 
undistinguished: ±t-\/|skp, which is exactly the disper- 
sion relation of a single layer for a given spin state32S!, 



For A 7^ 0, /i 7^ 0, we have ± 



(MTi|sfe|) for 



each layer, which is the spectrum for the s-wave pair- 
in#2l. 

If we take into account the hopping term between 
planes, in the absence of superconductivity, we obtain 
±|i?BL|- In particular, at ^ = 0, the four energy 
bands along three directions in the first Brillouin zone 
for A = /X = can be seen in Fig. [2ja, which is the 
same plot shown for the unbiased graphene bilayer irPSl. 
The eight distinct energy bands in the normal state are 
shown in Fig. [2]b. (It should be noticed that, for this 
particular choice of parameters, the chemical potential 
is sitting right at the bottom of the upper band.) As 
expected, the system is gapped in the superconducting 
state and the four energy bands at /i = are shown in 
Fig[2]c. Finally, the eight energy bands for nonzero gap 
and chemical potential are shown in Fig. [2jd. 

The graphene dispersion relation has six Dirac points 
at the corners of the first Brillouin zone; however, only 
two of them are non-equivalent. The continuum limit 
of our model Hamiltonian is obtained expanding Eq. ([s]) 
in the vicinity of the Dirac points K = —^'nj'^aex £^nd 
K' = 47r/3ae:r, 
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where a ^ K, K' and Aa is obtained replacing tsk by 
—vy {kx — iky) and —vp [k^ + iky) in Eq. ([s]) for K and 
K' respectively, with h = 1 and vp = \/3to/2. 

The partition function in the complex time represen- 
tation is written as 



Z = 



1 




, (13) 



where is the vacuum functional, 13 ~ l/ksT (fc^ is 
the Boltzmann constant) and 



rCL 
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with -f/)^ 



E 



d^k 
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"a,<T,; (^) ' ^I,<T,; (^) representing the spino- 
rial fields appearing in the continuum limit of the tight- 
biding graphene Hamiltonian density. 

Integrating over the fermion fields, we get that the par- 
tition function is proportional to 



idri^a. 
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an— icOnl+Aa is a function of the Matsubara 



where A'^ 

frequencies for fermions, w„ 
8x8 unity matrix, with ks 
of simplicity. 

Finally, redefining the coupling g 



2n + l)nT, and 1 is the 
1 hereafter for the sake 



\/N, the "effective 
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potential" per bilayer for each Dirac point will be 
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This would be the leading order in a 1 jN expansion and 
would be the exact result for TV — cx). 

In the next section, we analyze the conditions for 
the appearance of superconductivity at zero temperature 
provided by our mean-field effective potential. 



A. The superconducting instabilities at T = 

We shall study the minima of the effective potential. 
The occurrence of superconductivity corresponds to the 
existence of nonzero solutions of the order parameter, 
which minimizes the effective potential. 

Taking the derivative of Vefi with respect to the order 
parameter and summing over the Matsubara frequencies, 
we obtain 



K'ff(r) - A* 
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The nonzero solutions for |A| are given equalizing to 
zero the expression within brackets in the Eq. ( 17) above, 
which provides a self-consistent gap equation. In partic- 
ular, at zero temperature, we get 



Kff(O) = A* 



2 

A 2 
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(27r) y^A^Ti 



(19) 



Introducing a large momentum cutoff A/wp, we can 
integrate Eq. (19 1 over fc, 
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where a = iirvi., a,b — ±1, 



(20) 



(21) 



and €ab = Cab (A 0). 

In particular, for /x = t_L = 0, the nonzero solutions for 
the superconducting gap ar^^ll 



A = 



~2 



A2 



1 

A2 



(22) 



for A > a/ A, what establishes quantum critical point 
for the onset of superconductivity in the system at the 
critical coupling Ac = a/ A. 

Evidencing A in ( 20 ) , it can be reexpressed as 



off 




, (23) 



where A' = A/ Ac and the tilde indicates that the quan- 
tity is divided by A. Since all the nonzero solutions for 
the gap are given by the expression between the curly 
brackets above, and given that it does not depend on 
any experimental data, our results are suitable to de- 
scribe any planar Dirac fermion system with a hopping 
between adjacent sheets, assuming that A is the single 
free parameter of our model. Therefore, in the following 
we present our numerical results for A in terms of the 
parameter A. 

Notice that the cutoff is always provided by the lat- 
tice in condensed matter systems. Indeed, we have 
A = 2'KhvY/a as an upper bound for the energy cutoff, 
and since a is the smallest distance scale, A becomes a 
natural high-energy cutoff. In fact, this frequently hap- 
pens in condensed matter. A familiar example in the 
case of conventional, phonon mediated superconductiv- 
ity, is the Debye frequency (energy) a natural cutoff that 
emerges in BCS theory. Moreover, we also constrain our- 
selves to positive values of the chemical potential up to 
/i/A — 0.9, given the half bandwidth of A. 

The case for = Q and finite /i with different val- 
ues of interaction coupling has been exhaustively inves- 
tigated by some of us'^^. The plots of the superconduct- 
ing gap as a function of /i for A/ Ac = 0.8 are shown in 
Fig. |3] Starting at /i/A = 0, the system is in the normal 
state, since A < Ac. As /i/A increases, Aq/A displays 
a dome-shaped plot: the system asymptotically becomes 
superconducting up to a maximum value at an optimal 
chemical potential and decreases as /i/A increases even 
further. 

Notice that the system is not quantum critical but, in 
fact, the curve vanishes exponentially as /i — > 0, hence, 
superconductivity persists down to /i = 0. 

Our results are consistent with^^, which also obtains 
a dome-shaped plot of A for relativistic interacting par- 
ticles, as can be seen in Fig. 1 of their paper (choice of 
parameters I, referred as the weak- coupling case). 

An interesting result is obtained as we increase the 
value of t± , as can bee seen in the inset of Fig. [3j For 
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Figure 3: The superconducting gap as a function of the chemi- 
cal potential for t± = 0, 0.3 and 0.5. The inset shows the same 
plot for a smaller range of the chemical potential. A/Ac ~ 0.8 
and all the other quantities are given in unities of A. 



ter becomes nonzero, as determined by the gap equation. 
Therefore, Tc calculated in this section may be regarded 
as a mean-field upper bound critical temperature for the 
KT transition, which sets the actual temperature for the 
appearance of superconductivity in the N ^ oo limit. 
W e st art considering the gap equation provided by the 



Eq. ( 17|. Making the change of variables x 
get 



we 
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where 



= 0, (24) 



(25) 



with £,ab{x) given by Eq. (21), replacing for x. From 
Eq. ( [24| , we calculate the superconducting critical tem- 
perature Tc making A = at T = Tc in the above ex- 
pression. 



small values of the chemical potential, as the out of the 
plane hopping between layers increases, we see that the 
superconducting gap also increases for the same value of 
the chemical potential. Indeed, even for /i = 0, for which 
there is no superconducting gap when t± = 0, given that 
A < Ac, there is a nonzero value of A for t±/A = 0.3 or 
0.5, indicating that the system is in the superconducting 
state. Therefore, the hopping between layers favors the 
appearance of superconductivity. 

As shall be seen in the next section, since the energy 
gap and the superconducting critical temperature tend 
to be proportional quantities, this result also shows that 
the superconducting critical temperature increases as the 
hopping between layers increases for small values of the 
chemical potential. 



B. Superconducting phase at finite temperatures 

In this section we calculate the superconducting phase 
diagram for finite temperatures. A priori, the nonzero 
solutions for A are supposed to hold only in the N ^ oo 
limit at a finite temperature, because, otherwise, they are 
ruled out by the Coleman-Mermin-Wagner-Hohcnberg 
theorerrPni_ This limit corresponds to a physical situa- 
tion where the three-dimensionality of the system is ex- 
plicitly taken into account. For finite values of N and 
T ^ 0, there is an underlying Berezinskii-Kosterlitz- 
Thouless (BKT) transitiori2il, below which phase coher- 
ence is found for a nonzero A. The actual superconduct- 
ing transition occurs at Tbkt i£ ^c. However, it can be 

shown that Tbkt ^^^^ T^. This clearly indicates that, 
in spite of the fact that we may have a nonzero supercon- 
ducting gap at T = Tc, only in a really three-dimensional 
system we will have phase coherence developing at the 
same temperature that the modulus of the order parame- 
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Figure 4: The superconducting critical temperature as a func- 
tion of the chemical potential for t± = 0, 0.3 and 0.5. The 
inset shows the same plot for a smaller range of the chemical 
potential. A/ Ac = 0.8 and all the other quantities are given 
in unities of A. 



In Section |II A| we have found a dome-shaped super- 



conducting gap as a function of the chemical potential 
for A < 1 and several values of t^. Since the energy 
gap and the superconducting critical temperature Tc are 
proportional, we expect to find a dome-shaped plot for 
Tc as a function of as well. In fact, as can be seen 
in Fig. [4] our numerical results for the superconducting 
critical temperature presents the characteristic dome ex- 
perimentally observed in several compounds, like 1111 
pnictides and cuprate superconductors. 

A dome- like structure of the superconducting phase for 
two-dimensional Dirac fermions has been previously ob- 



tained irpS', where the superconducting critical temper- 
ature also presents a dome at intermediate filling frac- 
tions, surrounded by the normal phase for fillings close 
to unity or zero, which is consistent to our results. Also, a 



6 



dome for Tc as a function of hole concentrations has been 
previously obtained in a brief letter by some of u^^ for 
a relativistic version of the spin-fermion HamiltoniarPSl, 
used to describe the Cu-0 planes in the cuprates. Those 
results and the phase diagram presently calculated sug- 
gest that Dirac fermions may play a relevant role in the 
description of cuprates and iron pnictides. 

Indeed, it has been shown that Dirac points appear in 
the intersection of the nodes of the d-wave superconduct- 
ing gap and the 2D-Fermi surface in the high-Tc cuprate 
superconductors and the low-energy excitations will cor- 
respond exclusively to these pointpSl. Also, it has been 
experimentally found that the iron pnictide^^^'^Hl also 
present electronic excitations whose properties are gov- 
erned by the Dirac equation. Theoretical results also sup- 
port the existence of Dirac electrons in the pnictides'^^'^. 

As also suggested in the former section, the inset of 
Fig. [3] shows that, for a small value of /i, Tc increases as 
the i_L parameter is increased, indicating that the hop- 
ping between layers favors the appearance of supercon- 
ductivity in the system. 

As shall be seen in the next section, the same feature 
is observed as we take into account the first-neighbors 
out of the plane hopping between adjacent layers, which 
is the case for graphite. 



III. GRAPHITE 

In this section, we calculate the superconducting phase 
diagram of many coupled graphene layers for a finite 
chemical potential. To simplify the problem, we consider 
only the minimal model where only the electron tunnel- 
ing amplitudes between the nearest sites in the plane t 
and out of the plane t±_ are regarded. The same approach 
was employed irP3, as briefly explained below. 

Consider a Bernal-stacked graphene bilayer described 
by the following Hamiltonian in the vicinity of each non- 
equivalent Dirac poinlP^, 



k.cr 



where the above 4x4 matrix B]^ is given by 



6w = 



(26) 



(27) 



the vector a = {(7x,<^y) is written in terms of the well- 
known Pauli matrices, the matrix B12 is 



B12 = Sji = 



tj_ 




(28) 



introducing the operator 



k,cr 



the Hamiltonian becomes 



(29) 
(30) 



k,(T 



where 



vpli-a B12 

B21 WFk-CT B12 

B21 vp'k- a 



\ 



(31) 
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Introducing the momentum fc^ in the z direction, it is 
possible to re-express the Hamiltonian for graphite taking 
the first-neighbors hopping between adjacent layers in the 
momentum representation, which is wri tten in terms of 
a 4 X 4 matrix similar to Sk in Eq. ( 27 1^, 



k.fc,, 



where 



2?k„ 



vp^-S 2 B12 cos kzd 
2 B21 cos kzd UFk-CT 



(32) 



(33) 



and d is the distance between layers. 

For this minimal model, the dispersion relation is given 
by 



^Gr = ±y\vF k|2 + (t_L COS k,d) ± t±_ COS k^d (34) 

and, for kzd — 7r/2, we recover the Dirac-type dispersion 
found in graphene. 

Taking into account the attractive interaction form- 
ing Cooper pairs within each graphene layer, as seen in 
Eq. ([2|, and introducing the operator 



k.cr 



(35) 



where "^k o- ; given by Eq. ii» , the model Hamiltonian 
which describes the superconducting graphite in a mean- 
field approximation becomes 



k.CT 



(36) 



where 



and $L = 



'k,o-j- "k,cr,i 



('/'lc,aa''/'U,2)' with 
j = 1,2 denotes the layer index. 

The model Hamiltonian for graphite is assumed to be 
described as an infinite number of graphene layers cou- 
pled by the hopping between adjacent sheets. Therefore, 



Ai A12 
A21 A2 A12 
A21 Ai 



(37) 
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Figure 5: The superconducting critical temperature as a func- 
tion of the chemical potential for (a) t±_ — 0.3 and (b) t± — 0.5 
for both graphene bilayer (solid line) and graphite (dotted 
line). The inset in each pane shows the same plot for a smaller 
range of the chemical potential. A/Ac = 0.8 and all the other 
quantities are given in unities of A. 



with Ai = Ai and A\2 — -^21 given by the Eqs. ([s]) and 
([9]) respectively. 

Accordingly, it is possible to re-express the Hamilto- 
nian for the superconducting graphite in terms of an 8 x 
8 matrix, which is similar to A^ in Eq. 



k.fcj 



where 



k.fc. 



Ai 2A12 COS kzd 
2A21 cos k^d A2 



and the dispersion is given by the 8 eigenvalues 



(38) 



(39) 



(40) 



with Eg,: given by Eq. (34) 



Therefore, the self-consistent equation for the super- 
conducting gap becomes 



2 
A 



2 ^ 



d 

1 

E+{Kk, 



d^k 
J2nf 



tanh 



where the four values of £'_|_(k, kz), labeled by the index 
j in the above expression, are given by Eq,- in Eq. ( 34 1 , 
in analogy to the discussions in the previous sections. 
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E+{k,kz) 



, (41) 



We calculate the critical temperature from Eq. (41) 
and our results are compared with Tc obtained for 
graphene bilayer, from the former section. Our results 
are shown in Fig. [5] and we see that there is not an en- 
hancement of Tf. for every range of chemical potential for 
a given value of t±. However, given tj_, we always find 
that the critical temperature for graphite is bigger than 
the Tc obtained for graphene bilayer for small values of 
the chemical potential, what demonstrates that the first 
neighbors hopping between adjacent sheets favors the su- 
perconductivity in the system. 



IV. CONCLUSIONS 

In conclusion, in the present paper we have derived the 
effective potential for a stack of graphene layers with a 
hopping between adjacent sheets and an on-site attrac- 
tive interaction between electrons in a mean-field approx- 
imation. 

For a single layer or two adjacent coupled layers of 
graphene, a remarkable result was obtained for the su- 
perconducting critical temperature as a function of the 
chemical potential: it displays a dome-shaped curve, as 
experimentally observed in several compounds, like 1111 
pnictides and cuprate superconductors. This result sug- 
gests that Dirac fermions may play a relevant role in the 
description of cuprates and iron pnictides, which shall be 
object of further investigation. Indeed, a dome- like struc- 
ture of the superconducting phase is in agreement with 
previous resul ts for strongly interacting two-dimensional 
Dirac fermion^^SEH^ pointed out in'^^, our results can 
also be experimentally realized with ultracold atoms in a 
two-dimensional optical square lattice. 

Finally, considering a minimal model for graphite, tak- 
ing into account only the tunneling amplitudes between 
the nearest sites in the plane and out of the plan^Sl^ 
we have compared the superconducting critical temper- 
ature for graphite and graphene bilayer. We have seen 
that the Tc calculated for graphite is bigger than the 
one for graphene bilayer for a small value of /i, what 
might explain why intrinsic superconductivity is observed 
in HOPG. 
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